fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)

mcmc=goMCMC(mdl,data,smp,wrm)

mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')


y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')


m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')
}

ex9

explanatory variable obs.x also has noise  
x~N(x0.sx)  
y~N(b0+b1*x0,s)  

ex8-0.stan

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m1[i],s);
  }
}

ex9.stan

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x10;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
  real<lower=0> sx;
  vector[N] x0;
}  
model{
  for(i in 1:N){
    x[i]~normal(x0[i],sx);
    y[i]~normal(b0+b1*x0[i],s);
  }
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x0[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] x1;
  vector[N1] y1;
  for(i in 1:N1){
    x1[i]=normal_rng(x10[i],sx);
    m1[i]=b0+b1*x10[i];
    y1[i]=normal_rng(m1[i],s);
  }
}
n=20
x0=runif(n,0,20)
x=rnorm(n,x0,2)
y=rnorm(n,x0*2+5,2)
qplot(x,y)

n1=10

#explanatory variable do not has noise
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-0.stan') 
fn(mdl,data)
## Initial log joint probability = -47080 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
##       21      -39.0024     0.0023489   0.000547583           1           1       77    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.2 seconds.
##  variable estimate
##    lp__     -39.00
##    b0         8.41
##    b1         1.68
##    s          4.26
##    m0[1]      9.31
##    m0[2]     33.85
##    m0[3]     24.79
##    m0[4]     38.10
##    m0[5]      7.71
##    m0[6]     25.79
##    m0[7]      5.55
##    m0[8]     27.47
##    m0[9]     41.37
##    m0[10]    36.98
##    m0[11]     4.55
##    m0[12]    11.87
##    m0[13]    33.88
##    m0[14]    19.90
##    m0[15]    22.36
##    m0[16]    34.77
##    m0[17]    18.98
##    m0[18]    39.96
##    m0[19]    29.58
##    m0[20]    19.55
##    y0[1]     12.18
##    y0[2]     37.30
##    y0[3]     21.40
##    y0[4]     39.16
##    y0[5]      9.97
##    y0[6]     28.36
##    y0[7]     14.26
##    y0[8]     27.45
##    y0[9]     43.77
##    y0[10]    40.46
##    y0[11]    -0.78
##    y0[12]    10.67
##    y0[13]    32.46
##    y0[14]    23.24
##    y0[15]    22.81
##    y0[16]    33.88
##    y0[17]    20.17
##    y0[18]    36.20
##    y0[19]    32.55
##    y0[20]    19.58
##    m1[1]      4.55
##    m1[2]      8.65
##    m1[3]     12.74
##    m1[4]     16.83
##    m1[5]     20.92
##    m1[6]     25.01
##    m1[7]     29.10
##    m1[8]     33.19
##    m1[9]     37.28
##    m1[10]    41.37
##    y1[1]     -8.70
##    y1[2]      8.55
##    y1[3]      9.77
##    y1[4]     22.88
##    y1[5]     20.64
##    y1[6]     29.07
##    y1[7]     24.51
##    y1[8]     30.62
##    y1[9]     36.82
##    y1[10]    44.40
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.4 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -39.19 -38.84 1.37 1.13 -41.88 -37.72 1.02      388      541
##    b0       8.51   8.52 1.99 1.91   5.28  11.72 1.02      201      210
##    b1       1.67   1.67 0.17 0.16   1.40   1.95 1.02      304      420
##    s        4.83   4.69 0.89 0.79   3.65   6.43 1.01     1100     1088
##    m0[1]    9.40   9.42 1.92 1.84   6.32  12.50 1.02      201      212
##    m0[2]   33.81  33.84 1.41 1.34  31.44  36.07 1.01     2420     1450
##    m0[3]   24.80  24.83 1.10 1.09  22.99  26.53 1.01      513      698
##    m0[4]   38.04  38.07 1.70 1.61  35.24  40.78 1.01     1618     1362
##    m0[5]    7.81   7.83 2.05 1.97   4.45  11.12 1.02      202      212
##    m0[6]   25.79  25.80 1.10 1.09  23.97  27.51 1.01      658      864
##    m0[7]    5.66   5.68 2.23 2.17   1.96   9.30 1.02      204      220
##    m0[8]   27.47  27.49 1.12 1.12  25.59  29.21 1.01     1008     1276
##    m0[9]   41.29  41.33 1.97 1.86  38.00  44.45 1.01     1131     1325
##    m0[10]  36.93  36.95 1.62 1.53  34.27  39.54 1.01     1806     1327
##    m0[11]   4.68   4.69 2.32 2.26   0.80   8.45 1.02      204      222
##    m0[12]  11.95  11.96 1.71 1.67   9.17  14.74 1.02      203      233
##    m0[13]  33.84  33.86 1.41 1.34  31.46  36.10 1.01     2414     1450
##    m0[14]  19.94  19.95 1.21 1.14  17.90  21.91 1.01      261      283
##    m0[15]  22.39  22.41 1.13 1.10  20.51  24.20 1.01      329      366
##    m0[16]  34.72  34.75 1.47 1.37  32.27  37.08 1.01     2254     1466
##    m0[17]  19.03  19.02 1.25 1.20  16.87  21.07 1.01      247      254
##    m0[18]  39.89  39.92 1.85 1.75  36.79  42.85 1.01     1304     1344
##    m0[19]  29.56  29.58 1.19 1.18  27.57  31.41 1.01     1852     1528
##    m0[20]  19.59  19.60 1.23 1.17  17.50  21.59 1.01      255      257
##    y0[1]    9.43   9.42 5.13 5.01   0.90  17.89 1.00      917     1560
##    y0[2]   33.79  33.86 5.11 4.83  25.15  41.94 1.00     1778     1778
##    y0[3]   24.79  24.89 4.93 4.98  16.85  32.80 1.00     1385     1464
##    y0[4]   38.14  38.16 5.33 4.96  29.21  46.74 1.00     1868     1814
##    y0[5]    7.74   7.64 5.18 5.16  -0.44  16.19 1.00      854     1592
##    y0[6]   25.83  25.94 5.11 5.06  17.70  34.01 1.00     2015     1888
##    y0[7]    5.76   5.68 5.49 5.24  -3.47  14.79 1.00      893     1354
##    y0[8]   27.28  27.26 5.22 4.83  18.83  35.69 1.00     1786     1965
##    y0[9]   41.28  41.26 5.29 5.15  32.55  49.91 1.00     1794     1968
##    y0[10]  36.91  36.74 5.15 5.05  28.63  45.51 1.00     1941     1845
##    y0[11]   4.74   4.61 5.45 5.17  -3.93  13.83 1.01      674     1579
##    y0[12]  12.07  12.12 5.28 5.09   3.29  20.76 1.00     1192     1727
##    y0[13]  33.81  33.72 5.31 5.15  25.16  42.30 1.00     1906     1818
##    y0[14]  19.89  19.82 5.00 4.67  11.77  28.39 1.00     1783     1973
##    y0[15]  22.24  22.14 5.03 4.86  14.01  30.36 1.00     2022     1892
##    y0[16]  34.70  34.76 5.10 4.73  26.33  42.68 1.00     2133     1938
##    y0[17]  18.94  19.09 4.89 4.67  10.73  26.71 1.00     1503     1482
##    y0[18]  39.96  39.94 5.13 4.93  31.48  48.06 1.00     1918     2058
##    y0[19]  29.57  29.61 5.00 4.81  21.52  37.88 1.00     1942     1703
##    y0[20]  19.62  19.55 5.11 4.95  11.33  27.94 1.00     1947     1780
##    m1[1]    4.68   4.69 2.32 2.26   0.80   8.45 1.02      204      222
##    m1[2]    8.74   8.76 1.97 1.88   5.55  11.92 1.02      201      210
##    m1[3]   12.81  12.81 1.65 1.61  10.13  15.50 1.02      205      236
##    m1[4]   16.88  16.87 1.37 1.33  14.58  19.12 1.02      224      244
##    m1[5]   20.95  20.97 1.17 1.12  19.00  22.85 1.01      283      306
##    m1[6]   25.02  25.03 1.10 1.08  23.20  26.74 1.01      544      706
##    m1[7]   29.08  29.10 1.17 1.17  27.09  30.90 1.01     1590     1524
##    m1[8]   33.15  33.18 1.37 1.31  30.83  35.33 1.01     2521     1513
##    m1[9]   37.22  37.24 1.64 1.54  34.53  39.86 1.01     1753     1346
##    m1[10]  41.29  41.33 1.97 1.86  38.00  44.45 1.01     1131     1325
##    y1[1]    4.78   4.60 5.54 5.19  -4.41  14.23 1.00      942     1084
##    y1[2]    8.84   8.86 5.40 5.16   0.13  17.50 1.00      874     1663
##    y1[3]   12.71  12.75 5.23 5.15   4.42  21.28 1.00      849     1331
##    y1[4]   16.71  16.73 5.15 4.91   8.21  24.97 1.00     1559     1695
##    y1[5]   21.01  20.96 5.08 4.81  12.87  29.54 1.00     1883     1950
##    y1[6]   25.03  24.99 4.95 4.69  16.87  33.21 1.00     1852     1712
##    y1[7]   29.11  29.12 5.10 4.67  20.58  37.21 1.00     2023     1872
##    y1[8]   33.06  33.02 5.20 4.99  24.44  41.33 1.00     2048     1945
##    y1[9]   37.05  37.05 5.11 5.05  28.91  45.23 1.00     1907     1786
##    y1[10]  41.43  41.27 5.42 5.21  32.61  50.18 1.00     2118     1801

#explanatory variable has noise
x10=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x10=x10)

mdl=cmdstan_model('./ex9.stan')
mcmc=goMCMC(mdl,data,wrm=500,smp=1000)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 1 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 2 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 3 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 4 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 2 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 3 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 4 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 1 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 3 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 4 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 1 finished in 0.4 seconds.
## Chain 3 finished in 0.4 seconds.
## Chain 4 finished in 0.4 seconds.
## Chain 2 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 2 finished in 0.5 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.4 seconds.
## Total execution time: 0.6 seconds.
seeMCMC(mcmc,'[m,x]')
##  variable   mean median    sd   mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -45.00 -48.12 13.28 11.60 -62.44 -19.59 1.04       82       61
##    b0       8.47   8.40  1.92  1.72   5.29  11.66 1.01      566      968
##    b1       1.67   1.68  0.17  0.14   1.41   1.96 1.01      407      937
##    s        3.32   3.50  1.69  1.84   0.49   6.01 1.08       41       38
##    sx       1.86   1.90  1.11  1.23   0.27   3.54 1.04       73      107
##    x0[1]    2.11   1.94  1.70  1.99  -0.21   4.78 1.03      105     1193
##    x0[2]   18.52  18.18  2.97  3.87  14.85  23.10 1.04       86      375
##    x0[3]    9.00   9.09  1.27  1.18   6.89  10.88 1.01      353     1752
##    x0[4]   18.57  18.39  1.52  1.24  16.53  21.21 1.01      286      545
##    x0[5]   -1.69  -1.38  1.76  1.58  -4.75   0.52 1.02      222      609
##    x0[6]   10.40  10.41  1.18  0.83   8.41  12.33 1.01     2270     1360
##    x0[7]   -0.97  -1.17  1.52  1.28  -3.08   1.40 1.02      227      953
##    x0[8]   10.33  10.45  1.39  1.41   8.03  12.39 1.02      172      977
##    x0[9]   19.37  19.38  1.36  1.03  17.24  21.60 1.01      521      682
##    x0[10]  16.66  16.71  1.29  0.96  14.67  18.70 1.01      701     1375
##    x0[11]  -2.08  -2.10  1.42  1.10  -4.38   0.11 1.01      410     1037
##    x0[12]   3.46   3.36  1.63  1.94   1.17   6.00 1.03      132     1493
##    x0[13]  14.70  14.72  1.28  1.07  12.66  16.73 1.01      428     1525
##    x0[14]   5.36   5.47  1.66  1.85   2.62   7.60 1.03      123     1060
##    x0[15]   8.62   8.59  1.13  0.94   6.77  10.50 1.00      767     1193
##    x0[16]  16.14  15.96  1.32  0.98  14.24  18.48 1.01      601      491
##    x0[17]   3.89   4.10  2.20  2.82   0.29   6.71 1.04      100      717
##    x0[18]  18.76  18.71  1.35  0.94  16.69  21.03 1.01     1083      856
##    x0[19]  11.94  12.03  1.35  1.19   9.81  13.91 1.02      245     2039
##    x0[20]   6.39   6.43  1.15  0.83   4.44   8.22 1.02     1467     1650
##    m0[1]   12.05  12.25  2.90  3.49   7.27  16.15 1.03      121      984
##    m0[2]   39.32  39.06  4.65  6.40  32.61  45.86 1.05       74      712
##    m0[3]   23.53  23.47  2.09  2.07  20.31  26.91 1.01      316     2501
##    m0[4]   39.43  39.64  2.30  2.24  35.51  42.81 1.01      343     1837
##    m0[5]    5.75   5.46  2.70  2.85   1.96  10.26 1.02      188     1365
##    m0[6]   25.86  25.83  1.93  1.51  22.84  28.99 1.02     3717     1838
##    m0[7]    6.92   7.15  2.50  2.31   2.51  10.75 1.01      325     1142
##    m0[8]   25.75  25.69  2.31  2.48  22.24  29.52 1.02      184     2505
##    m0[9]   40.77  40.65  2.25  1.86  37.06  44.47 1.02     1536     1814
##    m0[10]  36.28  36.14  2.14  1.87  33.00  39.82 1.01     1141     2220
##    m0[11]   5.10   5.27  2.32  1.88   1.09   8.82 1.02     1055     1477
##    m0[12]  14.29  14.54  2.78  3.11   9.67  18.24 1.03      137     1275
##    m0[13]  33.00  32.87  2.09  1.91  29.77  36.53 1.01      635     2381
##    m0[14]  17.49  17.52  2.64  3.22  13.70  21.67 1.03      119     2546
##    m0[15]  22.89  22.97  1.90  1.66  19.82  25.91 1.01     1286     2022
##    m0[16]  35.40  35.48  2.08  1.82  31.86  38.69 1.00     1025     2344
##    m0[17]  15.05  15.11  3.54  4.71  10.05  20.39 1.04       90     1443
##    m0[18]  39.76  39.84  2.18  1.77  36.11  43.33 1.01     2392     2274
##    m0[19]  28.41  28.30  2.17  2.04  25.15  32.00 1.02      266     2527
##    m0[20]  19.19  19.08  1.87  1.56  16.26  22.34 1.01     2294     2455
##    y0[1]   12.11  12.95  4.64  4.10   3.66  18.34 1.01      416     2034
##    y0[2]   39.39  40.19  5.82  6.57  28.84  46.79 1.03      124     1128
##    y0[3]   23.50  23.03  4.19  3.30  17.17  31.03 1.01     2515     3194
##    y0[4]   39.43  40.00  4.38  3.50  31.74  46.15 1.02     1601     1294
##    y0[5]    5.74   4.94  4.62  3.78  -0.94  13.80 1.01     1190     2883
##    y0[6]   25.91  25.91  4.24  3.18  18.86  33.27 1.02     4117     2972
##    y0[7]    7.02   7.52  4.45  3.56  -0.78  13.98 1.02     1609     2660
##    y0[8]   25.74  25.11  4.30  3.49  19.45  33.44 1.01      920     2042
##    y0[9]   40.88  40.62  4.33  3.42  33.92  48.28 1.02     3250     2678
##    y0[10]  36.32  36.03  4.18  3.36  29.60  43.56 1.02     3885     2613
##    y0[11]   5.13   5.37  4.35  3.29  -2.35  12.34 1.03     2857     2069
##    y0[12]  14.29  15.17  4.62  3.74   5.67  20.57 1.01      404     1851
##    y0[13]  32.95  32.62  4.25  3.37  26.31  40.33 1.02     2656     3015
##    y0[14]  17.48  16.61  4.57  3.74  11.18  25.96 1.01      672     1851
##    y0[15]  22.84  23.19  4.24  3.29  15.36  29.63 1.02     3672     3128
##    y0[16]  35.40  35.70  4.27  3.37  27.91  42.15 1.02     3333     3434
##    y0[17]  15.01  14.10  5.12  5.03   8.32  24.31 1.02      191     1746
##    y0[18]  39.67  39.73  4.33  3.42  32.37  46.61 1.03     3156     2625
##    y0[19]  28.46  27.89  4.35  3.41  21.87  36.18 1.02     1576     2671
##    y0[20]  19.23  18.99  4.17  3.18  12.54  26.39 1.03     3679     2730
##    m1[1]    4.65   4.56  2.25  1.98   0.91   8.34 1.00      527      959
##    m1[2]    8.71   8.63  1.90  1.71   5.58  11.86 1.01      569      968
##    m1[3]   12.78  12.72  1.58  1.38  10.20  15.38 1.01      653      901
##    m1[4]   16.84  16.77  1.31  1.15  14.71  19.04 1.01      833     1150
##    m1[5]   20.90  20.82  1.13  0.99  19.10  22.79 1.01     1098     1314
##    m1[6]   24.97  24.92  1.06  0.95  23.26  26.74 1.01     1364     1424
##    m1[7]   29.03  29.01  1.15  1.05  27.21  30.96 1.01     1205     1454
##    m1[8]   33.10  33.13  1.35  1.21  30.90  35.26 1.02      937     1087
##    m1[9]   37.16  37.26  1.63  1.45  34.45  39.76 1.02      662      902
##    m1[10]  41.22  41.36  1.96  1.70  37.95  44.33 1.02      542      880
##    x1[1]   -2.27  -2.26  2.16  1.40  -5.75   1.26 1.02     4003     2466
##    x1[2]    0.13   0.15  2.12  1.42  -3.36   3.72 1.02     3776     1972
##    x1[3]    2.62   2.58  2.14  1.45  -0.84   6.24 1.02     4080     1881
##    x1[4]    4.96   4.98  2.17  1.36   1.35   8.63 1.02     3985     2895
##    x1[5]    7.42   7.42  2.21  1.44   3.82  11.01 1.02     3914     1734
##    x1[6]    9.84   9.85  2.15  1.39   6.37  13.41 1.02     3945     1714
##    x1[7]   12.27  12.29  2.15  1.44   8.68  15.85 1.02     4218     1823
##    x1[8]   14.76  14.72  2.15  1.44  11.36  18.41 1.01     3977     2317
##    x1[9]   17.14  17.18  2.23  1.41  13.57  20.77 1.02     4113     1912
##    x1[10]  19.60  19.59  2.22  1.47  15.83  23.22 1.02     3867     2322
##    y1[1]    4.63   4.42  4.40  3.64  -2.24  12.08 1.02     1642     2065
##    y1[2]    8.68   8.46  4.08  3.40   2.08  15.48 1.02     1730     2227
##    y1[3]   12.81  12.71  4.11  3.33   6.13  19.53 1.01     1907     2417
##    y1[4]   16.85  16.74  3.94  3.14  10.42  23.25 1.02     3393     2750
##    y1[5]   20.94  20.86  3.88  3.06  14.30  27.23 1.02     3606     2506
##    y1[6]   24.91  24.88  3.78  2.96  18.47  31.23 1.02     3973     3161
##    y1[7]   29.00  29.07  3.87  3.04  22.40  35.30 1.03     3454     1617
##    y1[8]   33.09  33.05  4.04  3.16  26.48  39.72 1.02     2508     2743
##    y1[9]   37.21  37.29  4.04  3.29  30.53  43.93 1.02     2510     2734
##    y1[10]  41.11  41.30  4.22  3.51  34.04  47.98 1.02     1593     2332

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
x1=mcmc$draws('x1')
smx1=summary(x1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=smx1$median,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

ex10

outlier
y~cauchy(b0+b1*x,s)

ex10.stan

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~cauchy(b0+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=cauchy_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=cauchy_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,9)
y=rnorm(n,x*2+5,1)
x[1]=3
y[1]=25
qplot(x,y)

n1=10

x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-0.stan') 
fn(mdl,data)
## Initial log joint probability = -1673.04 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       19      -32.8956     0.0011713   0.000402703      0.8809      0.8809       29    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__     -32.90
##    b0         7.11
##    b1         1.67
##    s          3.14
##    m0[1]     12.13
##    m0[2]     20.03
##    m0[3]     12.76
##    m0[4]      9.63
##    m0[5]     10.58
##    m0[6]     16.27
##    m0[7]     10.69
##    m0[8]     19.16
##    m0[9]     17.03
##    m0[10]     8.85
##    m0[11]    18.34
##    m0[12]    19.89
##    m0[13]    13.73
##    m0[14]    16.61
##    m0[15]    13.81
##    m0[16]    16.51
##    m0[17]     9.57
##    m0[18]    15.16
##    m0[19]    22.00
##    m0[20]    21.08
##    y0[1]     10.47
##    y0[2]     18.10
##    y0[3]     14.27
##    y0[4]     12.60
##    y0[5]     11.86
##    y0[6]     16.84
##    y0[7]      8.54
##    y0[8]     20.21
##    y0[9]     13.24
##    y0[10]     8.07
##    y0[11]    20.48
##    y0[12]    18.28
##    y0[13]    10.81
##    y0[14]    16.93
##    y0[15]    16.76
##    y0[16]    13.20
##    y0[17]     9.40
##    y0[18]    16.45
##    y0[19]    24.08
##    y0[20]    18.91
##    m1[1]      8.85
##    m1[2]     10.31
##    m1[3]     11.77
##    m1[4]     13.23
##    m1[5]     14.69
##    m1[6]     16.15
##    m1[7]     17.61
##    m1[8]     19.08
##    m1[9]     20.54
##    m1[10]    22.00
##    y1[1]      6.84
##    y1[2]     11.54
##    y1[3]      7.90
##    y1[4]     11.21
##    y1[5]     15.30
##    y1[6]     14.52
##    y1[7]     17.47
##    y1[8]     19.48
##    y1[9]     18.87
##    y1[10]    23.52
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -33.29 -32.98 1.22 1.02 -35.54 -31.95 1.01      662      878
##    b0       7.23   7.22 1.66 1.72   4.50   9.97 1.02      448      411
##    b1       1.65   1.66 0.32 0.32   1.15   2.18 1.01      515      543
##    s        3.53   3.46 0.64 0.60   2.65   4.69 1.00     1099     1340
##    m0[1]   12.18  12.19 0.94 0.96  10.58  13.72 1.01      595      796
##    m0[2]   19.99  19.99 1.25 1.17  17.92  22.04 1.00     1317     1450
##    m0[3]   12.81  12.81 0.88 0.90  11.33  14.25 1.01      677      829
##    m0[4]    9.72   9.72 1.26 1.29   7.57  11.81 1.01      472      435
##    m0[5]   10.65  10.66 1.13 1.14   8.75  12.52 1.01      498      536
##    m0[6]   16.27  16.26 0.83 0.80  14.94  17.65 1.00     2333     1397
##    m0[7]   10.77  10.77 1.11 1.13   8.89  12.60 1.01      503      555
##    m0[8]   19.12  19.12 1.13 1.07  17.27  21.00 1.00     1570     1428
##    m0[9]   17.02  17.02 0.89 0.86  15.58  18.47 1.00     2427     1308
##    m0[10]   8.94   8.94 1.38 1.41   6.64  11.25 1.02      459      429
##    m0[11]  18.32  18.31 1.02 0.98  16.64  19.98 1.00     1991     1421
##    m0[12]  19.85  19.85 1.23 1.16  17.83  21.86 1.00     1350     1464
##    m0[13]  13.76  13.76 0.82 0.83  12.39  15.10 1.00      884     1138
##    m0[14]  16.61  16.60 0.85 0.81  15.24  18.02 1.00     2401     1401
##    m0[15]  13.84  13.84 0.81 0.82  12.47  15.17 1.00      907     1151
##    m0[16]  16.51  16.51 0.85 0.80  15.15  17.91 1.00     2385     1422
##    m0[17]   9.65   9.65 1.27 1.30   7.49  11.76 1.01      471      433
##    m0[18]  15.17  15.18 0.79 0.77  13.85  16.47 1.00     1592     1384
##    m0[19]  21.92  21.92 1.56 1.50  19.31  24.48 1.01     1009     1407
##    m0[20]  21.02  21.01 1.41 1.34  18.65  23.29 1.00     1122     1471
##    y0[1]   12.32  12.25 3.71 3.55   6.33  18.43 1.00     1994     2000
##    y0[2]   20.04  20.05 3.86 3.63  13.69  26.40 1.00     1817     1805
##    y0[3]   12.90  12.88 3.64 3.51   6.95  18.60 1.00     1672     1930
##    y0[4]    9.69   9.76 3.82 3.69   3.35  15.73 1.00     1556     1829
##    y0[5]   10.56  10.60 3.79 3.69   4.16  16.68 1.00     1518     1705
##    y0[6]   16.17  16.20 3.72 3.47  10.21  22.22 1.00     2090     1877
##    y0[7]   10.66  10.67 3.67 3.63   4.87  16.67 1.00     1552     1872
##    y0[8]   19.23  19.23 3.69 3.59  13.21  25.24 1.00     2002     2042
##    y0[9]   17.00  17.09 3.64 3.45  11.08  22.92 1.00     2116     1926
##    y0[10]   8.82   8.91 3.79 3.50   2.66  15.04 1.00     1390     1916
##    y0[11]  18.20  18.24 3.76 3.70  12.05  24.34 1.00     2107     1969
##    y0[12]  19.63  19.74 3.82 3.74  13.26  25.84 1.00     2066     1848
##    y0[13]  13.81  13.87 3.70 3.74   7.74  19.93 1.00     1947     1864
##    y0[14]  16.65  16.76 3.74 3.74  10.39  22.66 1.00     1892     1794
##    y0[15]  13.71  13.74 3.81 3.62   7.50  19.94 1.00     1891     2080
##    y0[16]  16.69  16.66 3.68 3.58  10.55  22.83 1.00     1967     1964
##    y0[17]   9.56   9.57 3.88 3.64   3.14  15.95 1.00     1283     1313
##    y0[18]  15.22  15.29 3.71 3.56   8.87  21.38 1.00     2147     1954
##    y0[19]  21.85  21.93 3.99 3.65  15.27  28.31 1.00     1728     1651
##    y0[20]  21.04  21.08 3.91 3.66  14.55  27.33 1.00     1524     1848
##    m1[1]    8.94   8.94 1.38 1.41   6.64  11.25 1.02      459      429
##    m1[2]   10.39  10.39 1.17 1.18   8.41  12.32 1.01      489      507
##    m1[3]   11.83  11.85 0.98 0.99  10.17  13.42 1.01      563      692
##    m1[4]   13.27  13.27 0.85 0.86  11.85  14.65 1.01      760     1053
##    m1[5]   14.71  14.72 0.79 0.79  13.37  16.00 1.00     1296     1307
##    m1[6]   16.16  16.15 0.82 0.80  14.82  17.51 1.00     2297     1439
##    m1[7]   17.60  17.59 0.94 0.91  16.04  19.15 1.00     2381     1243
##    m1[8]   19.04  19.05 1.12 1.05  17.20  20.90 1.00     1604     1428
##    m1[9]   20.48  20.49 1.33 1.25  18.31  22.64 1.00     1214     1449
##    m1[10]  21.92  21.92 1.56 1.50  19.31  24.48 1.01     1009     1407
##    y1[1]    8.93   8.92 3.83 3.80   2.65  15.39 1.00     1547     1729
##    y1[2]   10.31  10.34 3.72 3.69   4.27  16.32 1.00     1613     1785
##    y1[3]   11.83  11.86 3.59 3.37   5.90  17.78 1.00     1616     1738
##    y1[4]   13.25  13.29 3.71 3.60   7.07  19.51 1.00     1884     1954
##    y1[5]   14.65  14.71 3.76 3.69   8.34  20.73 1.00     1719     1914
##    y1[6]   16.13  16.02 3.69 3.61  10.08  21.93 1.00     1899     2013
##    y1[7]   17.62  17.55 3.72 3.39  11.73  23.78 1.00     1987     2054
##    y1[8]   19.19  19.22 3.79 3.73  12.76  25.17 1.00     1944     1934
##    y1[9]   20.37  20.38 3.93 3.76  13.96  26.95 1.00     1842     1853
##    y1[10]  21.94  21.94 3.96 3.90  15.52  28.56 1.00     1918     1620

mdl=cmdstan_model('./ex10.stan') 
fn(mdl,data)
## Initial log joint probability = -118.839 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       21      -15.0002   7.45559e-05   0.000492157      0.8905      0.8905       35    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__     -15.00
##    b0         5.69
##    b1         1.84
##    s          0.75
##    m0[1]     11.20
##    m0[2]     19.89
##    m0[3]     11.90
##    m0[4]      8.46
##    m0[5]      9.50
##    m0[6]     15.75
##    m0[7]      9.62
##    m0[8]     18.92
##    m0[9]     16.59
##    m0[10]     7.59
##    m0[11]    18.03
##    m0[12]    19.73
##    m0[13]    12.96
##    m0[14]    16.12
##    m0[15]    13.05
##    m0[16]    16.01
##    m0[17]     8.39
##    m0[18]    14.53
##    m0[19]    22.04
##    m0[20]    21.03
##    y0[1]     11.85
##    y0[2]     20.25
##    y0[3]     10.37
##    y0[4]      7.82
##    y0[5]      9.65
##    y0[6]     16.38
##    y0[7]      9.46
##    y0[8]     18.58
##    y0[9]     12.31
##    y0[10]    10.16
##    y0[11]   -44.34
##    y0[12]    19.90
##    y0[13]    13.45
##    y0[14]    15.79
##    y0[15]    15.72
##    y0[16]    16.37
##    y0[17]     7.54
##    y0[18]    15.08
##    y0[19]    16.34
##    y0[20]    20.49
##    m1[1]      7.59
##    m1[2]      9.20
##    m1[3]     10.81
##    m1[4]     12.41
##    m1[5]     14.02
##    m1[6]     15.62
##    m1[7]     17.23
##    m1[8]     18.83
##    m1[9]     20.44
##    m1[10]    22.04
##    y1[1]      8.11
##    y1[2]      8.26
##    y1[3]     11.13
##    y1[4]     13.07
##    y1[5]     12.20
##    y1[6]     13.64
##    y1[7]     15.91
##    y1[8]     18.57
##    y1[9]     20.32
##    y1[10]    23.13
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median     sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -16.89 -16.60   1.29 1.08 -19.41 -15.45 1.00      579      688
##    b0       5.52   5.57   0.66 0.65   4.40   6.54 1.02      318      422
##    b1       1.85   1.85   0.11 0.10   1.68   2.03 1.01      377      515
##    s        0.87   0.84   0.26 0.24   0.51   1.36 1.00      854      677
##    m0[1]   11.08  11.11   0.42 0.41  10.35  11.73 1.01      416      693
##    m0[2]   19.85  19.85   0.41 0.40  19.17  20.50 1.00     2206     1597
##    m0[3]   11.79  11.81   0.39 0.39  11.10  12.41 1.01      461      776
##    m0[4]    8.31   8.35   0.53 0.53   7.39   9.12 1.02      335      550
##    m0[5]    9.36   9.40   0.48 0.49   8.51  10.11 1.01      355      606
##    m0[6]   15.68  15.69   0.33 0.35  15.13  16.20 1.00     1369     1417
##    m0[7]    9.49   9.52   0.48 0.48   8.65  10.23 1.01      358      606
##    m0[8]   18.87  18.87   0.38 0.38  18.25  19.48 1.00     2257     1508
##    m0[9]   16.52  16.52   0.34 0.36  15.96  17.05 1.00     1801     1471
##    m0[10]   7.44   7.49   0.57 0.56   6.47   8.32 1.02      327      511
##    m0[11]  17.97  17.98   0.36 0.37  17.38  18.53 1.00     2219     1515
##    m0[12]  19.69  19.69   0.41 0.40  19.02  20.34 1.00     2220     1546
##    m0[13]  12.86  12.88   0.36 0.36  12.22  13.43 1.01      583      842
##    m0[14]  16.05  16.06   0.33 0.35  15.51  16.57 1.00     1595     1471
##    m0[15]  12.94  12.97   0.36 0.36  12.32  13.51 1.01      597      828
##    m0[16]  15.94  15.95   0.33 0.35  15.39  16.46 1.00     1538     1481
##    m0[17]   8.24   8.28   0.53 0.53   7.32   9.06 1.02      334      550
##    m0[18]  14.44  14.46   0.34 0.35  13.87  14.97 1.01      863     1070
##    m0[19]  22.02  22.02   0.50 0.47  21.21  22.82 1.00     1504     1415
##    m0[20]  21.01  21.00   0.45 0.44  20.25  21.72 1.00     1934     1411
##    y0[1]   12.07  11.13  25.14 1.38   5.12  17.02 1.00     1918     1997
##    y0[2]   19.63  19.87  25.01 1.37  14.32  25.78 1.00     1903     1925
##    y0[3]   12.81  11.87  43.76 1.41   6.18  17.82 1.00     1790     1930
##    y0[4]    8.26   8.36  47.67 1.43   3.28  14.11 1.00     1706     1968
##    y0[5]    8.60   9.36  18.96 1.45   3.45  15.00 1.00     1693     1973
##    y0[6]   16.74  15.70  23.48 1.38  11.27  21.19 1.00     1946     2055
##    y0[7]    9.37   9.42  26.47 1.34   3.43  14.33 1.00     1674     1932
##    y0[8]   19.44  18.90  18.09 1.32  13.58  24.09 1.00     1922     1956
##    y0[9]   16.93  16.51  55.12 1.42  11.13  21.63 1.00     2021     1843
##    y0[10]   6.93   7.45  26.45 1.41   2.11  12.69 1.00     1684     1915
##    y0[11]  19.36  17.98  45.64 1.39  12.64  23.34 1.00     2039     1959
##    y0[12]  17.93  19.72  37.46 1.36  14.16  24.59 1.00     1888     1931
##    y0[13]  13.27  12.86  15.69 1.33   8.12  17.88 1.00     2102     2039
##    y0[14]  16.43  15.99  43.79 1.30  10.74  21.12 1.00     1852     1832
##    y0[15]   9.27  12.99 214.56 1.37   7.57  19.16 1.00     2095     1971
##    y0[16]  16.62  15.94  39.23 1.30  11.36  20.89 1.00     1949     1973
##    y0[17]   8.67   8.32  20.13 1.41   4.10  14.63 1.00     1693     2027
##    y0[18]  18.90  14.44 163.39 1.32   9.33  20.33 1.00     1905     1850
##    y0[19]  20.86  22.06  41.15 1.44  17.02  27.37 1.00     2052     2017
##    y0[20]  21.34  21.01  29.23 1.36  16.17  26.20 1.00     2016     1845
##    m1[1]    7.44   7.49   0.57 0.56   6.47   8.32 1.02      327      511
##    m1[2]    9.06   9.10   0.49 0.50   8.19   9.82 1.01      348      585
##    m1[3]   10.68  10.72   0.43 0.43   9.93  11.35 1.01      398      703
##    m1[4]   12.30  12.33   0.38 0.38  11.64  12.90 1.01      510      832
##    m1[5]   13.92  13.95   0.34 0.35  13.34  14.46 1.01      751      961
##    m1[6]   15.54  15.55   0.33 0.35  14.99  16.07 1.00     1295     1352
##    m1[7]   17.16  17.17   0.34 0.36  16.61  17.71 1.00     2081     1539
##    m1[8]   18.78  18.79   0.38 0.38  18.16  19.39 1.00     2256     1508
##    m1[9]   20.40  20.40   0.43 0.41  19.70  21.08 1.00     2139     1455
##    m1[10]  22.02  22.02   0.50 0.47  21.21  22.82 1.00     1504     1415
##    y1[1]    5.08   7.49  93.65 1.46   1.28  13.51 1.00     1901     1944
##    y1[2]    9.41   9.08  30.62 1.45   2.96  13.94 1.00     1737     2056
##    y1[3]   11.42  10.74  34.16 1.40   4.86  16.42 1.00     1984     1766
##    y1[4]    9.79  12.28  75.51 1.51   6.10  18.30 1.00     1895     2010
##    y1[5]   11.54  13.90 160.53 1.38   8.19  19.70 1.00     2033     2040
##    y1[6]   16.05  15.58  42.81 1.34  10.17  21.27 1.00     1874     1918
##    y1[7]   16.13  17.18  49.36 1.31  11.42  23.08 1.00     1873     1914
##    y1[8]   20.54  18.75  67.97 1.40  12.61  24.85 1.00     1854     1953
##    y1[9]   20.92  20.44  50.03 1.40  14.60  26.21 1.00     1841     1934
##    y1[10]  22.58  22.00  23.62 1.38  15.86  27.32 1.00     1904     1971

ex11

censored  

y i=1-N, y~N(m,s)  
  acutualy        ya i=1-Na
  lower censored  yl i=1-Nl, y<L, P(y<L)=cdf(L-m /s)
  upper censored  yu i=1-Nu, y>U, P(y>U)=ccdf(U-m /s)

cdf(z) cumulative normal density function, P((-Infinity,z],z~N(0,1))
ccdf(z) complementary CDF, P([z,Infinity),z~N(0,1))

P(y | x,m,s)=P(ya i=1-Na)* P(yl i=1-Nl)* P(yu i=1-Nu)

ex11.stan

data{
  int N;
  vector[N] x;
  vector[N] y;
  real L;
  int Nl;
  vector[Nl] xl;
  real U;
  int Nu;
  vector[Nu] xu;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0+b1*x,s);
  for(i in 1:Nl)
    target+=normal_lcdf(L | b0+b1*xl[i],s);
  for(i in 1:Nu)
    target+=normal_lccdf(U | b0+b1*xu[i],s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m1[i],s);
  }
}
n0=20
x=runif(n0,0,9)
y=rnorm(n0,10+3*x,4)
L=15
y[y<L]=L
nl=sum(y==L)
U=30
y[y>U]=U
nu=sum(y==U)
n=n0-nl-nu
qplot(x,y)

xy0=tibble(x=x,y=y)
xya=filter(xy0, y>L & y<U)
xyl=filter(xy0, y==L)
xyu=filter(xy0, y==U)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=xya$x,y=xya$y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan') 

mle=mdl$optimize(data=data)
## Initial log joint probability = -5182.13 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       21      -9.04902   0.000449008   1.46199e-05      0.9943      0.9943       44    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__      -9.05
##    b0        12.41
##    b1         2.19
##    s          1.66
##    m0[1]     19.49
##    m0[2]     20.64
##    m0[3]     19.32
##    m0[4]     16.70
##    m0[5]     17.44
##    m0[6]     16.77
##    m0[7]     18.42
##    m0[8]     26.00
##    m0[9]     16.65
##    y0[1]     20.39
##    y0[2]     19.66
##    y0[3]     19.25
##    y0[4]     16.92
##    y0[5]     16.74
##    y0[6]     16.10
##    y0[7]     18.49
##    y0[8]     24.81
##    y0[9]     16.58
##    m1[1]     12.73
##    m1[2]     14.85
##    m1[3]     16.96
##    m1[4]     19.07
##    m1[5]     21.19
##    m1[6]     23.30
##    m1[7]     25.42
##    m1[8]     27.53
##    m1[9]     29.64
##    m1[10]    31.76
##    y1[1]     11.41
##    y1[2]     15.89
##    y1[3]     14.11
##    y1[4]     17.12
##    y1[5]     20.59
##    y1[6]     22.58
##    y1[7]     25.16
##    y1[8]     28.41
##    y1[9]     31.86
##    y1[10]    31.81
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,'m')
##  variable   mean median   sd  mad     q5   q95 rhat ess_bulk ess_tail
##    lp__   -10.40  -9.98 1.61 1.27 -13.67 -8.72 1.01      364      333
##    b0      12.65  12.64 2.16 1.82   9.25 16.04 1.01      311      246
##    b1       2.12   2.15 0.65 0.54   1.05  3.14 1.01      338      307
##    s        2.32   2.12 0.86 0.64   1.38  3.91 1.01      770      603
##    m0[1]   19.51  19.53 0.87 0.72  18.12 20.80 1.00     1318     1050
##    m0[2]   20.62  20.63 0.98 0.81  19.04 22.06 1.00     1334      886
##    m0[3]   19.34  19.36 0.87 0.72  17.96 20.65 1.00     1229     1025
##    m0[4]   16.81  16.80 1.11 0.94  15.07 18.56 1.00      392      334
##    m0[5]   17.52  17.52 0.99 0.84  15.99 19.08 1.00      473      476
##    m0[6]   16.87  16.87 1.10 0.92  15.16 18.61 1.00      398      349
##    m0[7]   18.47  18.48 0.88 0.75  17.07 19.86 1.00      742      762
##    m0[8]   25.81  25.92 2.25 1.85  22.17 29.18 1.00      462      518
##    m0[9]   16.75  16.75 1.12 0.94  15.01 18.53 1.00      388      333
##    y0[1]   19.48  19.55 2.66 2.20  15.32 23.34 1.00     1739     1660
##    y0[2]   20.60  20.62 2.71 2.21  16.36 24.76 1.01     1706     1454
##    y0[3]   19.30  19.37 2.70 2.22  14.84 23.35 1.00     1941     1428
##    y0[4]   16.89  16.90 2.69 2.29  12.61 21.01 1.00     1235     1466
##    y0[5]   17.46  17.51 2.68 2.29  13.19 21.79 1.00     1508     1414
##    y0[6]   16.81  16.79 2.78 2.38  12.56 21.27 1.00     1089      998
##    y0[7]   18.48  18.44 2.63 2.20  14.31 22.92 1.00     1783     1552
##    y0[8]   25.81  25.86 3.28 2.83  20.56 30.98 1.00      826     1083
##    y0[9]   16.69  16.70 2.66 2.29  12.42 20.93 1.00     1329     1061
##    m1[1]   12.96  12.95 2.07 1.76   9.69 16.23 1.01      312      240
##    m1[2]   15.01  15.00 1.52 1.29  12.63 17.43 1.00      326      276
##    m1[3]   17.06  17.05 1.06 0.89  15.40 18.74 1.00      416      394
##    m1[4]   19.10  19.13 0.86 0.72  17.73 20.42 1.00     1089      944
##    m1[5]   21.15  21.18 1.07 0.89  19.41 22.76 1.00     1099      817
##    m1[6]   23.20  23.28 1.53 1.25  20.68 25.46 1.00      603      580
##    m1[7]   25.24  25.34 2.08 1.73  21.86 28.38 1.00      480      518
##    m1[8]   27.29  27.40 2.67 2.21  22.94 31.34 1.01      433      465
##    m1[9]   29.34  29.49 3.28 2.71  23.94 34.38 1.01      410      470
##    m1[10]  31.39  31.56 3.89 3.21  24.99 37.37 1.01      395      452
##    y1[1]   13.05  12.94 3.40 2.87   7.77 18.69 1.00      598      704
##    y1[2]   15.05  14.96 2.98 2.50  10.58 19.99 1.00      837     1087
##    y1[3]   17.13  17.08 2.67 2.41  12.97 21.52 1.00     1151     1168
##    y1[4]   19.12  19.10 2.67 2.21  14.98 23.49 1.00     1805     1581
##    y1[5]   21.06  21.12 2.75 2.34  16.55 25.36 1.00     1977     1586
##    y1[6]   23.23  23.32 2.95 2.57  18.47 27.65 1.00     1208     1152
##    y1[7]   25.21  25.20 3.20 2.80  20.20 30.24 1.00      939     1049
##    y1[8]   27.30  27.42 3.73 3.21  21.31 32.99 1.00      652      684
##    y1[9]   29.34  29.61 3.98 3.45  22.76 35.30 1.00      630      667
##    y1[10]  31.30  31.49 4.65 4.13  23.99 38.21 1.00      517      609

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

data=list(N=n,x=xya$x,y=xya$y,
          L=L,Nl=nl,xl=xyl$x,
          U=U,Nu=nu,xu=xyu$x,
          N1=n1,x1=x1)
mdl=cmdstan_model('./ex11.stan') 

mle=mdl$optimize(data=data)
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Initial log joint probability = -263.167 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
##       26      -15.7809   0.000214419   0.000130948           1           1       36    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -15.78
##    b0         9.07
##    b1         3.14
##    s          2.44
##    m0[1]     19.22
##    m0[2]     20.86
##    m0[3]     18.97
##    m0[4]     15.23
##    m0[5]     16.29
##    m0[6]     15.33
##    m0[7]     17.69
##    m0[8]     28.55
##    m0[9]     15.15
##    y0[1]     18.42
##    y0[2]     17.31
##    y0[3]     18.82
##    y0[4]     13.22
##    y0[5]     15.08
##    y0[6]     12.99
##    y0[7]     22.07
##    y0[8]     30.35
##    y0[9]     15.84
##    m1[1]      9.54
##    m1[2]     12.57
##    m1[3]     15.60
##    m1[4]     18.63
##    m1[5]     21.65
##    m1[6]     24.68
##    m1[7]     27.71
##    m1[8]     30.74
##    m1[9]     33.76
##    m1[10]    36.79
##    y1[1]     14.98
##    y1[2]     11.13
##    y1[3]     15.60
##    y1[4]     12.19
##    y1[5]     19.70
##    y1[6]     25.30
##    y1[7]     30.84
##    y1[8]     29.33
##    y1[9]     35.88
##    y1[10]    30.31
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
seeMCMC(mcmc,'m',ch=T)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -16.59 -16.16 1.40 1.15 -19.41 -15.07 1.00      491      630
##    b0       7.86   8.18 2.30 2.06   3.59  11.00 1.01      289      370
##    b1       3.52   3.42 0.64 0.56   2.68   4.82 1.01      311      430
##    s        3.36   3.08 1.14 0.88   2.04   5.52 1.01      617      572
##    m0[1]   19.25  19.27 1.03 0.93  17.54  20.92 1.00      962      861
##    m0[2]   21.09  21.04 1.09 0.92  19.38  22.94 1.00     1443     1206
##    m0[3]   18.97  18.99 1.03 0.93  17.26  20.59 1.00      881      854
##    m0[4]   14.76  14.92 1.31 1.19  12.48  16.62 1.01      365      509
##    m0[5]   15.95  16.07 1.19 1.08  13.89  17.66 1.01      424      529
##    m0[6]   14.88  15.03 1.29 1.18  12.62  16.71 1.01      369      517
##    m0[7]   17.54  17.60 1.07 0.97  15.74  19.13 1.00      586      690
##    m0[8]   29.72  29.43 2.18 1.84  26.81  34.01 1.01      582      618
##    m0[9]   14.68  14.83 1.32 1.21  12.39  16.53 1.01      362      509
##    y0[1]   19.18  19.27 3.69 3.29  13.16  25.05 1.00     1916     1522
##    y0[2]   21.16  21.13 3.68 3.25  15.32  27.19 1.00     1968     1716
##    y0[3]   18.99  18.97 3.64 3.22  13.19  24.72 1.00     1475     1505
##    y0[4]   14.86  15.08 3.78 3.45   8.46  20.55 1.00     1253     1572
##    y0[5]   16.01  15.95 3.76 3.21  10.11  21.92 1.00     1455     1412
##    y0[6]   14.78  14.90 3.69 3.30   8.81  20.52 1.00     1318     1281
##    y0[7]   17.59  17.54 3.83 3.44  11.74  23.42 1.00     1820     1722
##    y0[8]   29.82  29.53 4.07 3.60  23.85  37.02 1.00     1136     1101
##    y0[9]   14.65  14.79 3.85 3.34   8.38  20.64 1.00     1252     1008
##    m1[1]    8.38   8.69 2.21 1.98   4.24  11.40 1.01      291      370
##    m1[2]   11.78  11.99 1.69 1.52   8.53  14.14 1.01      307      383
##    m1[3]   15.18  15.32 1.26 1.16  12.98  16.98 1.01      382      525
##    m1[4]   18.58  18.61 1.04 0.93  16.86  20.18 1.00      786      853
##    m1[5]   21.98  21.91 1.15 0.98  20.26  24.02 1.00     1535     1048
##    m1[6]   25.38  25.21 1.53 1.30  23.25  28.27 1.00     1070      871
##    m1[7]   28.78  28.51 2.03 1.71  26.05  32.69 1.01      636      695
##    m1[8]   32.18  31.82 2.58 2.20  28.79  37.24 1.01      492      602
##    m1[9]   35.58  35.13 3.16 2.67  31.45  41.93 1.01      430      593
##    m1[10]  38.98  38.43 3.75 3.17  34.09  46.48 1.01      397      583
##    y1[1]    8.34   8.65 4.38 3.84   0.90  14.74 1.00      798     1102
##    y1[2]   11.83  11.93 3.89 3.43   5.19  17.86 1.00      909     1383
##    y1[3]   15.08  15.17 3.63 3.32   9.10  20.84 1.00     1421     1560
##    y1[4]   18.57  18.69 3.66 3.33  12.35  24.24 1.00     1486     1500
##    y1[5]   21.91  21.84 3.71 3.22  16.01  28.04 1.00     1918     1693
##    y1[6]   25.50  25.29 4.00 3.43  19.51  32.09 1.00     1634     1590
##    y1[7]   28.95  28.61 4.15 3.66  23.02  36.21 1.00     1186     1118
##    y1[8]   32.11  31.71 4.41 3.82  25.74  40.00 1.00     1206     1115
##    y1[9]   35.64  35.16 4.71 4.20  28.89  43.91 1.00     1018     1254
##    y1[10]  38.80  38.17 5.12 4.66  31.58  48.23 1.00      672      825

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

bimodal distribution
mixed normal distribution
data {
  int<lower=0> N;               // データの数
  real y[N];                    // 観測データ
}
parameters {
  real<lower=0, upper=1> theta; // 分布1と分布2の混合比率
  real mu1;                     // 分布1の平均
  real mu2;                     // 分布2の平均
  real<lower=0> sigma1;         // 分布1の標準偏差
  real<lower=0> sigma2;         // 分布2の標準偏差
}
model {
  for (n in 1:N) {
    target += log_mix(theta,
                      normal_lpdf(y[n] | mu1, sigma1),
                      normal_lpdf(y[n] | mu2, sigma2));
  }
}
y0=rnorm(100,0,1)
y1=rnorm(100,-5,1)
y2=rnorm(100,5,1)
y=sample(c(y0,y1,y2),100)
density(y) |> plot()

EM algorithm

library(mclust)

rst=Mclust(y)
summary(rst)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust E (univariate, equal variance) model with 3 components: 
## 
##  log-likelihood   n df  BIC  ICL
##            -246 100  6 -520 -520
## 
## Clustering table:
##  1  2  3 
## 29 33 38
rst$parameters
## $pro
## [1] 0.290 0.331 0.378
## 
## $mean
##      1      2      3 
## -4.977  0.234  4.842 
## 
## $variance
## $variance$modelName
## [1] "E"
## 
## $variance$d
## [1] 1
## 
## $variance$G
## [1] 3
## 
## $variance$sigmasq
## [1] 0.925
## 
## 
## $Vinv
## NULL
plot(rst)